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Abstract: The El Nino phenomenon, synonymously El Nino-Southern Oscillation (ENSO), is 
an anomalous climatic oscillation in the Equatorial Pacific that occurs once every 3-8 years. 
It affects the earth's climate on a global scale. Whether it is a cyclic or a sporadic event, 
or whether its apparently random behaviour can be explained by stochastic dynamics have 
remained matters of debate. Herein ENSO is viewed, unconventionally, as a two-dimensional 
dynamical system on a desktop. The main features of ENSO: irregularity, interannual vari- 
ability, and the asymmetry between El Nino and La Nina are captured, simply, comprehen- 
sibly, quickly and cheaply, in a nonlinear stochastic limit cycle paradigm. Its predictability 
for ENSO compares remarkably well with that of the best state-of-the-art complex models 
from the European and American meteorological centres. Additionally, for the first time, 
by analyzing subsurface Equatorial Pacific data since 1960, this model finds that long-term 
variations are not caused by ENSO itself, but by external sources. 
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El Nino is an anomalous climatic event in the Equatorial Pacific. On average, the western 
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Pacific sea-surface temperature (SST) is about 4°-5°C higher than that of the eastern Pacific, and 
the thermocline, or the sharp vertical temperature gradient separating the warm and the cold water, 
lies deeper in the western Pacific than in the Eastern Pacific. The higher SST of the western Pacific 
results in rising warm, moist air from the Indonesian coast, and the descent of cool, dry air towards 
the eastern south Pacific. These atmospheric movements are supplemented by the high altitude 
eastward jet stream and the westward trade winds along the Equatorial Pacific sea surface. During 
an El Nino, which peaks in the northern winter, the eastern Pacific gets warmer as the pool of warm 
water stretches much further eastward. The western Pacific thermocline becomes shallower and 
the westward trade winds at the sea surface become much weaker, and sometimes even reverse. 
These anomaly conditions in the Equatorial Pacific affect the earth's climate at a global scale. 

That ENSO is a result of ocean-atmosphere interaction has been clear since the time of 
Bjerknes 1 . It has been well-established that there are three key players in the ocean- atmosphere 
interaction that determine the ENSO dynamics at interannual time-scales: the SST, the atmospheric 
wind stress anomaly and the ocean's ability to transport energy 2 in the Equatorial Pacific. However, 
how the dynamics of these quantities relate to each other to perturb the normal climatic conditions 
to create an El Nino, and especially how the growth of an El Nino is controlled so that the normalcy 
is restored have remained matters of discussion. 3 A large number of models geared to address 
these issues — ranging from dynamical models inspired by the physics of ENSO 416 to stochastic 
ones 17 27 — have put forward different mechanisms at varying degrees of complexity: for some 
ENSO is inherently unstable 45 and the nonlinearities in ENSO dynamics eventually control the 
growth of an El Nino 6-11 , while some others have assumed that the normal climatic condition 
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is stable and that ENSO is excited by noise 1517 ' 20 27 . These models have also differed in their 
interpretation of the role of noise originating from external processes at a faster time scale, and 
the role of nonlinearity in ENSO properties. Over and above, they can capture different features 
of ENSO in their own ways; nevertheless, the ability to reproduce all the basic features of ENSO, 
while being able to simultaneously predict ENSO with good accuracy in a comprehensive manner 
is still lacking in the ENSO model world 3 . Remarkably, one such feature of ENSO not easily 
reproduced is the skewness of the ENSO indices, or the asymmetry between El Nino and La Nina 28 . 

Given this setting, herein ENSO is viewed, unconventionally, as a two-dimensional dynam- 
ical system on a desktop, described by the Nino3 index, the anomaly in the average SST of the 
region 5°S-5°N, 150°W-90°W, and the anomaly in Z20 or the average depth of the 20°C isotherm 
(as a proxy for the thermocline depth) of the region 5°S-5°N, 120°E-290°E, henceforth denoted 
by T and H respectively. Using the Nino3 and Z20 data averaged over month i, the (discrete) 
ENSO dynamical system indexed by (T i: Hi) is constructed. Based on the physics of ENSO, a 
phenomenological model, subject to fixed- amplitude Gaussian white noise, is conjectured to de- 
scribe this dynamical system. The model captures, simply, comprehensibly, quickly and cheaply, 
all the main features of ENSO: irregularity, interannual variability, and the asymmetry between El 
Nino and La Nina, in a nonlinear stochastic limit cycle paradigm. Its predictability for ENSO com- 
pares remarkably well with that of the best state-of-the-art complex models from the European and 
American meteorological centres. Additionally, for the first time, by analyzing subsurface Equa- 
torial Pacific (Z20) data since 1960 (roughly the time when regular record-keeping of subsurface 
Pacific temperature profile began), this model finds that decadal (and above) variations stem from 
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external sources, e.g. climate shifts. 

ENSO as a nonlinear stochastic limit cycle 

The idea behind viewing ENSO as a two-dimensional dynamical system is based on the so-called 
"recharge oscillator" 15 ' 29 ' 30 , which identifies the Equatorial eastern Pacific SST anomaly T E , the 
equatorial eastern and the western Pacific thermocline depth anomalies H E and H w , and the cen- 
tral Pacific zonal wind-stress anomaly r as the main oceanic and atmospheric quantities involved in 
the ENSO dynamics. The Equatorial eastern and the western Pacific thermocline depth anomalies 
can be combined to form two independent variables: (H E — H w ), representing the "thermocline 
anomaly tilt", and (H E + H w ), representing the anomaly in the amount of warm water volume 
present in the Equatorial Pacific (WWVA) 15,29-32 . Of these two, a positive anomaly in (H E — H w ) 
is strongly in phase with a positive eastern Pacific SST anomaly and a positive central Pacific 
zonal wind-stress anomaly. As envisaged by the recharge oscillator, before the onset of an El 
Nino, (meridional) Sverdrup transport of warm water towards the Equator gives rise to a positive 
anomaly in the WWVA. A positive WWVA gives rise to a positive T E . Subsequently, r responds 
positively to a positive T E , and increases (H E — H w ). These changes in the Pacific surface con- 
ditions then generate Rossby and Kelvin thermocline waves in the equatorial waveguide: these 
waves first make the western Pacific thermoclines shallower, and the return of the reflected waves 
at the western Pacific boundaries, further on, reduces the eastern Pacific SST anomaly. 

The recharge oscillator does not explain how the (meridional) Sverdrup transport of warm 
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water towards (or away from) the Equator is triggered. Nonetheless, it seems logical that a positive 
WWVA results in a positive T E and vice versa. A linear relationship between these two quantities 
was conjectured by Burgers et al. 15 , but that is in contradiction with the observation data: the 
magnitude of the positive T E due to a given positive magnitude of the WWVA during an El Nino 
is larger than the magnitude of T E due to a negative WWVA of the same magnitude during a La 
Nina 32 - 33 . This asymmetry between the El Nino and the La Nina (manifested equivalently via the 
Nino3 skewness 28 ), not fully understood at present, indicate that how the WWVA and T E interact 
is unclear: scenarios based on air-sea fluxes feedback on the SST, and the oceanic upwelling and 
vertical mixing mechanism have been suspected to contribute to it 34,35 . The fact however remains 
that the WWVA leads (H E — H w ), T E and r, which are strongly in phase with each other, by 
approximately 7 months on average, and consequently it should be considered as an important 
predictor for T E at ENSO time-scales. The natural variables in the (minimalistic) model for ENSO 
herein, therefore, are the Nino 3 monthly anomaly as a representative of T E , and the Z20 monthly 
anomaly, representing the WWVA (hereafter denoted by T and H respectively), forming a two- 
dimensional dynamical system. The lack of our understanding of how the WWVA and T E interact 
motivated this study of ENSO as a stochastic dynamical system, based on empirical data. 

The monthly Nino3 and Z20 anomalies in dimensionless units (rendered dimensionless by 
normalizing w.r.t. their r.m.s. magnitudes), hereafter denoted by (t, h), are shown in Fig. 1. 
A positive t (resp. h) means a positive Nino3 anomaly (resp. WWVA) and vice versa. In this 
notation, the aim of this study is to describe the ENSO dynamics stochastically as t i+1 = /(tj, hi) + 
6> hi+i — g(ti, hi) + rji, where f(x, y) and g(x, y) are two nonlinear functions of their arguments, 
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and £ and r] are mutually uncorrected fixed-amplitude Gaussian white noise. Visual inspection 
of the data in Fig. 1 suggests strong locking of the ENSO dynamics to its phase [defined by 
= — tan _1 (/i/t)], and therefore a cylindrical co-ordinate system r = \/t 2 + h 2 and 0, is more 
suitable. In these co-ordinates, without any loss of generality, the nonlinearities in the ENSO 
dynamics are then expressed by two coupled nonlinear differential equations as: dr/dt = r^f(r, fa) 
and d(f)/dt = uj(r, <fi), with the (r i: fa) -values obtained from the corresponding (t u hi) ones. 
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Figure 1: (a) Normalized monthly anomaly data for Nino3 (t) and Z20 (h) in the Equato- 
rial Pacific, describing ENSO for 1980-2004 (source: MET office objective analysis for the 
ENSEMBLES project, hereafter referred to as the "ENSEMBLES data"; publicly available 
at http://www.ecmwf.int/research/EU_projects/ENACT/ocean_analyses/index.html); (b) The same 
for the 1980-2002 data obtained from Australian Bureau of Meteorology (hereafter ABOM); pub- 
licly available at http://www.bom.gov.aU/bmrc/ocean/results/climocan.htm#subsurface%20anals]. 

Based on the resemblance of Fig. 1 to a limit cycle (i.e., a tendency to grow when the system 
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is close to r = 0, combined with a tendency to decay when the system grows far away form r = 0) 
strongly locked to its phase, the ansatz 7 (r, 0) = fi(<f>)[f2(<j>) —r] andcj(r, 0) = (yf 1 (0)y / ris made 
to describe ENSO. It is clear from Fig. 1 that /2(0) cannot simply be a constant, and for an expla- 
nation for uj(r, 0) oc ^Jr see the methods section. With this ansatz, the natural choice for /i(0), 
f 2 {4>) and gi(4>) is clearly a series expansion in increasing orders of sines and cosine harmonics as 
(A + Bi cos + B 2 sin + C 1 cos 20 + C 2 sin 20 + . . .). The parameters A, Bi, B 2 , . . . can then 
be estimated from a time-series (U, hi). This estimation is performed in the following manner: 
the observed time-series for N sequential months were denoted as (t 1: hi), (t 2 ,h 2 ), . . . ,(t N ,h N ). 
Then for i = 1, 2 . . . N, (t i: hi) is taken as the starting value and dr/dt and d(j)/dt are integrated 
forward from time % to time {% + 1) using the functional forms of 7(7", 0) and uj(r, 0). This yields 
the model theoretical values (t i+1 , h i+ i), as well as the noise = (t i+ i — t i+1 , h i+1 — h i+ i). The 
parameters A, Bi, B 2 , . . . can then be estimated by least-square optimization, i.e., by minimizing 
E = J^iLi 1 e i ' e i- N° te that this optimization process not only yields the parameter values, but 
also the characteristics of the noise e^, to be later used to realistically construct £ and 77. 

To fix the functional forms for /i(0), /2(0) and g(4>) that are to be used throughout the 
entire length of this study, ABOM data [Fig. 1(b)] were used for the least-square optimization. To 
capture the annual variations of t and h, any choice of /i(0), ^2(0) and g(4>) needed to include 
the fourth harmonics for the sine and cosine functions (on average, ENSO has a 4-year period); 
and simultaneously, due to the finiteness of the time-series, the number of parameters in 7(0) and 
u{4>) have been kept as low as possible. Eventually, the functional forms of /i(0), ^2(0) and g(4>) 
that used the smallest number of parameters while still leading to the smallest value of E for the 
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Figure 2: Phase-space manifold structure for the ENSO dynamics, given by dr/dt and d<p/dt, in 
terms of 7(7", <f>) and uj(r, 0).The parameters that generated these plots have been optimized on 
the 1980-2004 (t, h) ENSEMBLES data [Fig. 1(a)]. The corresponding phase-space manifold 
structure for the ABOM data (not shown here) of Fig. 1(b) is only marginally different from the 
above, as it should be; this confirms the stability of the optimization process to describe ENSO. 

Australian dataset were chosen. This procedure showed that there are 16 parameters needed to 
describe f\(4>), f2{4>) an d g{4>) [° r l{ r , 4>) an d u(r, </>)]. See the methods section for details. 

As it turns out, the ansatz about the functional forms for 7(7", 0) and u>(r, <p) are very well- 
chosen. One of its immediate outcomes is that the noise characteristics that emerged from Fig. 1 
can be fairly accurately described as Gaussian white. The rest of the outcomes: "stability" of this 
method, ENSO irregularity, variability, Nino3 skewness and predictability are discussed below. 
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Figure 3: (a): Probability density for finding an ENSO state at (t, h) in colour plot. The white 
closed curve is the average path (obtained by averaging the locations of all data points within 
10° ^-intervals for the 5, 000, 000 year run) of ENSO as a stochastic limit cycle. It implies that on 
average an initial ENSO state inside the white curve evolves in time to converge to it from inside, 
while an initial ENSO state outside the white curve evolves in time to converge to it from outside. 
The resemblance of Fig. 3(a) to Fig. 1(a) or (b) is clear. Note the relation of the white curve to 
the WWVA anomaly at the onset of an El Nino: a sharp positive anomaly in Equatorial WWVA 
precedes an El Nino, which is nicely captured by the hump in the white curve in the (t < 0, h > 0) 
quadrant, (b) Probability density P(s) of skewness s of the Nino3 index, calculated from 200, 000 
sets of 25 years each (obtained from the same 5, 000, 000 year run). The Nino3 skewness based 
on the 25 year UK MET office objective analysis data (1980-2004) appearing in Fig. 1(a) is 
approximately 0.85, while the average Nino3 skewness obtained from Fig. 3(b) is 0.81. The broad 
probability distribution of the Nino3 skewness exemplifies the extent of ENSO variability within 
the scope of this model. These plots use the same values for the parameters as in Fig. 2. 
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Stability of the optimization procedure 

The first issue for describing ENSO by means of these 16-parameter functional forms for 7(7", 0) 
and uj(r, 0) is the "stability" of the optimization procedure; namely that all the major ENSO fea- 
tures: irregularity, variability and Nino3 skewness are expected to be properly reproduced when 
the same minimization procedure is applied to a dataset reasonably similar to the ABOM one. The 
noise characteristics emerging from that new dataset are also expected to appear to be sufficiently 
close to Gaussian white. To this end the ENSEMBLES data [Fig. 1(a)] were used, and both expec- 
tations were fulfilled. Not only did the phase-space manifold structure of ENSO dynamics show 
clear signatures of growth and decay in the ENSO magnitude r [see Figs. 2(a) and (b)], but also 
when the corresponding parameters were used, in addition to Gaussian white noise (of strength 
obtained from the minimization procedure; see methods section), to simulate the two-dimensional 
stochastic dynamical system version of ENSO on a computer for 5, 000, 000 years, it produced the 
right 1980-2004 ENSO features. Two of these are shown in Fig. 3. In Fig. 3(a) appears the proba- 
bility density function for finding the ENSO state at (t, h) in colour plot. The white closed curve is 
the average path (obtained by averaging the locations of all data points within 10° ^-intervals for 
the 5, 000, 000 year run) of ENSO as a stochastic limit cycle. It implies that on average an initial 
ENSO state inside the white curve evolves in time to converge to it from inside, while an initial 
ENSO state outside the white curve evolves in time to converge to it from outside. The resem- 
blance of Fig. 3(a) to Fig. 1(a) or (b), which prompted the ansatz to model ENSO as a limit cycle 
in the first place, is clear. Also clear is the relation of the white curve to the WWVA anomaly at the 
onset of an El Nino: a sharp positive anomaly in Equatorial WWVA precedes an El Nino 32 ' 34 ' 35 , 
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which is nicely captured by the hump in the white curve in the (t < 0, h > 0) quadrant. Simulta- 
neously, Fig. 3(b) shows the probability density P(s) of skewness s of the Nino3 index, calculated 
from 200, 000 sets of 25 years each (obtained from the same 5, 000, 000 year run). The skewness 
of the Nino3 index based on the 25 year (1980-2004) ENSEMBLES data appearing in Fig. 1(a) is 
approximately 0.85, while the average Nino3 skewness obtained from Fig. 3(b) is 0.81. The broad 
probability distribution of the Nino3 skewness exemplifies the extent of ENSO variability within 
the scope of this model. 

Predictability skill 

The ultimate test of how robustly ENSO can be described by a two-dimensional stochastic dynam- 
ical system as above is to be able to predict ENSO with a reasonable accuracy. In order to address 
this issue, predictability for Nino3 over the period 1960-now was studied by running a 10,000 
different sequences of Gaussian white noise realizations in this model, using the ENSEMBLES 
data (available up to 1960, roughly when records-keeping of reliable regular subsurface Pacific 
temperature profile began). It is to be emphasized here that (a) The use of subsurface Equatorial 
Pacific data to analyze the relation between WWVA and ENSO has so far dated back to 1980. This 
study, therefore, is the first one to extend that to pre-1980; (b) This model's predictability for Nino3 
has been studied thoroughly by considering a wide variety of "training periods", i.e., the periods 
from which the data are used to optimize the parameters of the model. The results presented here 
are based on the training period 1980-2004 as described below; this choice is motivated as it also 
sheds light on the role of climate shifts on ENSO; and (c) Unless otherwise stated, in order to avoid 
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Figure 4: Comparison between the predicted average (red) and the observed (blue) values of nor- 
malized Nino3 monthly anomalies over the period 1960-now (and beyond): lead time 3 months 
(top) and 6 months (bottom). The yellow areas show the 2cr-spread in the 3- and 6-month lead time 
predictions (they also give a general idea about the magnitude of the 2cr-spread over the 1960-June 
2006 period), based upon runs of 10,000 realizations. See methods section for details. 

artificial skill, the target period has never been included in the training period. 

The predictability skill results (see methods section for the calculations and the convention 
for lead time) for the years 1960-now appear in Fig. 4. The underlying subtleties can be broken 
down to three separate periods, (i) 1980-2004: In order to make sure that the training period of the 
model did not include any information from the "target dataset", i.e., the dataset corresponding to 
the year for which Nino3 is to be predicted, the so-called "jackknifing" procedure was used. More 
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explicitly, the (t, h) normalized anomaly dataset was first formed out of the 1979-2004 Nino3 and 
Z20 monthly averages. Next, for the predictions of year n (1980 < n < 2004), the data for the 
years n — 1, n and n + 1 were removed from the (t, h) dataset of 1979-2004 to form a jackknived 
dataset, which served as the training period for the model, (ii) 1960-1979: Two separate normalized 
anomaly datasets (t, h) were formed, one out of the Nino3 and Z20 monthly averages of 1960- 
2004 data, and the other out of the 1980-2004 data. The target dataset was formed by truncating 
the former to 1960-1979, while the training period for the model was 1980-2004. The comparison 
between the corresponding predicted and observed t-values for the years 1960-1979 revealed a 
startling gap between the two (see Fig. 6 in Appendix D), as if the predicted magnitudes of Nino 3 
were uniformly shifted upwards by approximately 1°C compared to the observed ones over the 
period 1960-1976. Interestingly, when the Nino3 and Z20 anomaly values for 1960-2004 were 
first passed through a high-pass filter that removed the variations in these variables at time-scales 
> 10 years, and were subsequently normalized and truncated to 1960-1980 in order to form the 
target dataset, the uniform gap between the predicted and observed Nino3 values disappeared (this 
is the comparison shown in Fig. 4). Upon further reflection, it was understood that this uniform 
gap corresponds to the well-known climatological shift in 1976 (known as the "1976-shift" 42 ). 
Interpreted differently, the disappearance of the gap between predicted and filtered observation 
data indicates that variations in the Equatorial Pacific climatic conditions ENSO have two distinct 
components in it: climatological shifts (i.e., shifts in the mean background) that occur due to 
dynamics at decadal (or above) time-scales and variations that occur at ENSO time-scales. It is the 
latter kind of variations that are captured by the stochastic limit cycle picture. The issues related 
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Figure 5: (a) Predictability skill of this model for 1960-June 2006 (i.e., correlations between the 
red and blue curves of Fig. 4). See also Fig. 7 in Appendix D for skill as a function of months, (b) 
Nino3 monthly anomaly for 2006, based upon the climatology of 1980-2004; observation (blue); 
predicted average (red). The yellow area shows the 2<r-spread. 

to the origin of decadal variations in ENSO do remain matters of discussion 37 ^ 1 ; nevertheless, the 
conclusions of this analysis contradicts the view that climatological shifts of ENSO are a part of 
its own dynamics, (iii) 2005-now (and beyond): The climatological means and normalizations for 
the target dataset (t, h) as well as the training period for the model is 1980-2004. 

Finally, the predictability skill, i.e., Nino3 monthly anomaly correlations between observa- 
tion and average of Nino3 anomaly (as they appear in Fig. 4) for 1960-June 2006 are shown in Fig. 
5(a). The predictions until December 2006, as obtained by (iii) above appear in Fig. 5(b). 
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Skill comparison with operational models 



The comparison of the predictability skill of this model with the four highest-ranked operational 
models [System-2 (S2) 43 coupled atmosphere-ocean model from European meteorological centre 
ECMWF, and Constructed Analogue (CA) 44 ' 45 , Markov 46 and Climate Forecast System (CFS) 47 
from the U.S. centre for meteorology NCEP] for the period 1987-2001 are presented in Table I. 
It provides a fair yardstick for the predictability power of this model. Since different strategies 
were used in different models to generate ensembles and the number of members varied over the 
period 1987-2001, making this comparison has not been entirely trivial. Nevertheless, the period 
1987-2001 is actually selected as it is the common denominator of these operational models 48 . See 
also Fig. 8 in Appendix D for skill comparisons as a function of months. 
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0.76 


1987-2001 


Markov 
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0.83 
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0.90 


0.79 
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Table I: Comparison of the prediction skill of this dynamical systems model with the best opera- 
tional models [S2 from ECMWF, and CA, Markov and CFS from NCEP] for the period 1987-2001. 
The CA and Markov models are statistical, and the CFS and S2 models are coupled atmosphere- 
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ocean general circulation model (AO-GCM). Over this period, the two-dimensional dynamical sys- 
tem model used here performs slightly better than those of the NCEP statistical ones, and slightly 
worse than the AO-GCM of ECMWF. Source: van Oldenborgh et al, 2003 48 . 

It is also worth noting that a relatively recent work by Chen et al. 49 , using training period 
1980-2000 50 reported that large El Ninos could be predicted up to two years in advance. A pre- 
dictability comparison with that work before 1960 is not possible because of the lack of Z20 data. 
Nevertheless, when the training period 1980-2004 is used, the model discussed in this paper pro- 
duces very similar skills on target periods 1960-1975 and on 1976-1995 (Fig. 9 in Appendix D). 

Conclusion 

The idea behind viewing ENSO as a two-dimensional stochastic dynamical system in this work 
has been based on the recharge oscillator formulation. In this formulation, essentially the amount 
of warm water in the Equatorial Pacific controls ENSO. Although it is clear that a larger amount of 
warm water volume in the Equatorial Pacific results in a higher eastern Pacific SST and vice versa, 
how these two quantities interact is far from clear; a fact that this study draws its motivation from. 
It considers subsurface Pacific data roughly since record-keeping of reliable regular subsurface 
Pacific temperature profile began (1960), for the first time. This simple model firmly establishes 
that a two-dimensional dynamical system, and more specifically, a nonlinear stochastic limit cycle 
captures the essentials of ENSO to a great detail, as it (i) captures the well-known features of ENSO 
comprehensively, (ii) is able to predict ENSO with an accuracy well-comparable with the state-of- 
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the-art complex models from the European and American meteorological centres, and (iii) is able 
to separate the nature of decadal and multi-decadal variations from the interannual variability of 
ENSO. This development signifies the fact that oceanic upwelling, vertical layer mixing, and air- 
sea fluxes feedback mechanisms 34,35 should be paid more attention to in order to fully understand 
ENSO, as the interaction between the subsurface warm water volume and the SST takes place via 
these mechanisms. The expectation is that the parameters of this dynamical system are related to 
these mechanisms, and that remains a topic of future research. 

Appendix A 

Explanation for u(r, 0) oc y/r : This proportionality has been actually chosen (i) to avoid the 
ambiguity of defining u(r, 0) at r = 0, (ii) as well as for "regularization": should the system get 
stuck in a phase of uncontrolled growth, the \fr- term serves to move it quickly away from there. In 
general, u(r, 0) oc r n for any n > can serve for both, so there is some arbitrariness in the choice 
of n. Nevertheless, n — 1/2 was a choice motivated to keep the model as simple as possible: 
n = does not work for (i) and (ii), and n — 1 was found to contradict the observation data, and 
hence the choice n — 1/2 was made. It is also worthwhile to note that the results presented here 
are insensitive of the value of n within the range 0.3 < n < 0.8. 

Appendix B 

Least-square optimization and the functional forms of /i(0), /2(0) an d g{<t>) '• For given func- 
tional forms of fi((f>), ^2(0) and g(<f>) parametrized by n quantities (such as A, B\ etc.), the 
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least-square optimization is a minimization of the scalar function E defined on an n-dimensional 
manifold. In order to perform this minimization the amoeba method was used 36 . To make sure 
that the true minimum for E was reached, the minimization was started from a multitude of ini- 
tial values of the n parameters. A single minimization run on 25 years of monthly data takes 
about 15 minutes on a 1.8 GHz CPU. For the 1980-2002 (t, h) ABOM data initial trial minimiza- 
tion runs were set up with f x (0) = a + (6 01 cos + c i sin 0) + . . . + (b 04 cos 40 + c 4 sin 40), 
^2(0) = 1+&icos0+Ci sin 0, and g(0) = a 2 +(b 2 i cos0+c 2 i sin0) + . . .+(624 cos40+c 2 4 sin 40). 
It was found that increasing the number of parameters beyond 16 hardly reduced the minimum 
value of E. Once the number of parameters was thus fixed at 16, more trials were run with dif- 
ferent combinations of sine and cosine harmonics. Eventually, for the ABOM data, it was found 
that the forms of /i(0), /2(0) and g(4>) that lead to the minimum value of E are the following: 
/i(0) = a + &oicos0 + 6 2 cos 20 + c O2 sin20 + fe O 3COs3(0 - O 3) + & O 4Cos4(0 - O 3/2), 
^2(0) = 1 + b\ cos + ci sin 0, and g(0) = a 2 + 621 cos + c 2 i sin + b 2 \ cos 20 + c 2i sin 20 + 
6 23 cos 30 + 6 24 cos 40. These functional forms were used all throughout this study. The noise 
properties of obtained from these data were found to be Gaussian white to a good approximation. 

For later references, note here that when these functional forms were used on the 1980-2004 
ENSEMBLES data [see Fig. 1(a)], the noise characteristics again turned out Gaussian white to a 
good approximation, with HWHM 0.34 and 0.29 for t and h respectively. 



18 



Appendix C 



ENSO prediction as a function of lead months: The predictability studies for ENSO were per- 
formed in the usual way: the numerical values of the model parameters were obtained from the 
(t, h) data over the training period. These parameters were then used to time-evolve an initial 
state (t , h ) to (t i: hi) for % — 1, . . . , 12 using 10, 000 different sequences of Gaussian white noise 
realizations; this takes only about a minute on a 1.8 GHz CPU. The HWHM used for the Gaus- 
sian white noise were 0.34 for t and 0.29 for h, as found from the optimization procedure on the 
ENSEMBLES data. The predicted (U, hi) values constitute the prediction for "lead time % months". 
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Appendix D 

This appendix consists of additional figures to supplement the text. 




1960 1962 1964 1966 1968 1970 1972 1974 1976 1978 



Figure 6: Comparison between predicted (red) and unfiltered observation data (blue) for Nino3 
monthly anomaly over 1960-1979. The uniform gap of nearly 1°C between the two curves before 
1976 is due to the well-known 1976-shift. This gap disappears (Fig. 4 in the paper) when the 
observation data is passed through a high-pass filter that removes decadal (and above) variations 
in Nino3 and Z20. 
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Figure 7: Predictability skill of the dynamical system model as a function of target months for 



1960-2006 (i.e., between the red and blue curves of Fig. 4). 
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Figure 8: Predictability skill of the dynamical system model in the paper as a function of target 
months. This figure is meant for comparison with Fig. 2 of van Oldenborgh et al. (2005), which 
shows the ENSO predictability skill for the operational models of Table I. 
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Figure 9: Predictability skill of the dynamical systems model in the paper. This figure is meant 
for comparison with Fig. 2 of Chen et al. (2004). To generate this figure, the decadal (and above) 
variations have been removed from the 1960-1979 dataset, as discussed in the "Predictability skill" 
section of the paper. 
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